<!DOCTYPE html>
<html>
<head>
<title>时间序列与隐Markov链.md</title>
<meta http-equiv="Content-type" content="text/html;charset=UTF-8">

<style>
/* https://github.com/microsoft/vscode/blob/master/extensions/markdown-language-features/media/markdown.css */
/*---------------------------------------------------------------------------------------------
 *  Copyright (c) Microsoft Corporation. All rights reserved.
 *  Licensed under the MIT License. See License.txt in the project root for license information.
 *--------------------------------------------------------------------------------------------*/

body {
	font-family: var(--vscode-markdown-font-family, -apple-system, BlinkMacSystemFont, "Segoe WPC", "Segoe UI", "Ubuntu", "Droid Sans", sans-serif);
	font-size: var(--vscode-markdown-font-size, 14px);
	padding: 0 26px;
	line-height: var(--vscode-markdown-line-height, 22px);
	word-wrap: break-word;
}

#code-csp-warning {
	position: fixed;
	top: 0;
	right: 0;
	color: white;
	margin: 16px;
	text-align: center;
	font-size: 12px;
	font-family: sans-serif;
	background-color:#444444;
	cursor: pointer;
	padding: 6px;
	box-shadow: 1px 1px 1px rgba(0,0,0,.25);
}

#code-csp-warning:hover {
	text-decoration: none;
	background-color:#007acc;
	box-shadow: 2px 2px 2px rgba(0,0,0,.25);
}

body.scrollBeyondLastLine {
	margin-bottom: calc(100vh - 22px);
}

body.showEditorSelection .code-line {
	position: relative;
}

body.showEditorSelection .code-active-line:before,
body.showEditorSelection .code-line:hover:before {
	content: "";
	display: block;
	position: absolute;
	top: 0;
	left: -12px;
	height: 100%;
}

body.showEditorSelection li.code-active-line:before,
body.showEditorSelection li.code-line:hover:before {
	left: -30px;
}

.vscode-light.showEditorSelection .code-active-line:before {
	border-left: 3px solid rgba(0, 0, 0, 0.15);
}

.vscode-light.showEditorSelection .code-line:hover:before {
	border-left: 3px solid rgba(0, 0, 0, 0.40);
}

.vscode-light.showEditorSelection .code-line .code-line:hover:before {
	border-left: none;
}

.vscode-dark.showEditorSelection .code-active-line:before {
	border-left: 3px solid rgba(255, 255, 255, 0.4);
}

.vscode-dark.showEditorSelection .code-line:hover:before {
	border-left: 3px solid rgba(255, 255, 255, 0.60);
}

.vscode-dark.showEditorSelection .code-line .code-line:hover:before {
	border-left: none;
}

.vscode-high-contrast.showEditorSelection .code-active-line:before {
	border-left: 3px solid rgba(255, 160, 0, 0.7);
}

.vscode-high-contrast.showEditorSelection .code-line:hover:before {
	border-left: 3px solid rgba(255, 160, 0, 1);
}

.vscode-high-contrast.showEditorSelection .code-line .code-line:hover:before {
	border-left: none;
}

img {
	max-width: 100%;
	max-height: 100%;
}

a {
	text-decoration: none;
}

a:hover {
	text-decoration: underline;
}

a:focus,
input:focus,
select:focus,
textarea:focus {
	outline: 1px solid -webkit-focus-ring-color;
	outline-offset: -1px;
}

hr {
	border: 0;
	height: 2px;
	border-bottom: 2px solid;
}

h1 {
	padding-bottom: 0.3em;
	line-height: 1.2;
	border-bottom-width: 1px;
	border-bottom-style: solid;
}

h1, h2, h3 {
	font-weight: normal;
}

table {
	border-collapse: collapse;
}

table > thead > tr > th {
	text-align: left;
	border-bottom: 1px solid;
}

table > thead > tr > th,
table > thead > tr > td,
table > tbody > tr > th,
table > tbody > tr > td {
	padding: 5px 10px;
}

table > tbody > tr + tr > td {
	border-top: 1px solid;
}

blockquote {
	margin: 0 7px 0 5px;
	padding: 0 16px 0 10px;
	border-left-width: 5px;
	border-left-style: solid;
}

code {
	font-family: Menlo, Monaco, Consolas, "Droid Sans Mono", "Courier New", monospace, "Droid Sans Fallback";
	font-size: 1em;
	line-height: 1.357em;
}

body.wordWrap pre {
	white-space: pre-wrap;
}

pre:not(.hljs),
pre.hljs code > div {
	padding: 16px;
	border-radius: 3px;
	overflow: auto;
}

pre code {
	color: var(--vscode-editor-foreground);
	tab-size: 4;
}

/** Theming */

.vscode-light pre {
	background-color: rgba(220, 220, 220, 0.4);
}

.vscode-dark pre {
	background-color: rgba(10, 10, 10, 0.4);
}

.vscode-high-contrast pre {
	background-color: rgb(0, 0, 0);
}

.vscode-high-contrast h1 {
	border-color: rgb(0, 0, 0);
}

.vscode-light table > thead > tr > th {
	border-color: rgba(0, 0, 0, 0.69);
}

.vscode-dark table > thead > tr > th {
	border-color: rgba(255, 255, 255, 0.69);
}

.vscode-light h1,
.vscode-light hr,
.vscode-light table > tbody > tr + tr > td {
	border-color: rgba(0, 0, 0, 0.18);
}

.vscode-dark h1,
.vscode-dark hr,
.vscode-dark table > tbody > tr + tr > td {
	border-color: rgba(255, 255, 255, 0.18);
}

</style>
<link rel="stylesheet" href="file:///Users/william/Folders/mybooks/statistical-learning-abc/chapters/purple.css" type="text/css">
<style>
/* Tomorrow Theme */
/* http://jmblog.github.com/color-themes-for-google-code-highlightjs */
/* Original theme - https://github.com/chriskempson/tomorrow-theme */

/* Tomorrow Comment */
.hljs-comment,
.hljs-quote {
	color: #8e908c;
}

/* Tomorrow Red */
.hljs-variable,
.hljs-template-variable,
.hljs-tag,
.hljs-name,
.hljs-selector-id,
.hljs-selector-class,
.hljs-regexp,
.hljs-deletion {
	color: #c82829;
}

/* Tomorrow Orange */
.hljs-number,
.hljs-built_in,
.hljs-builtin-name,
.hljs-literal,
.hljs-type,
.hljs-params,
.hljs-meta,
.hljs-link {
	color: #f5871f;
}

/* Tomorrow Yellow */
.hljs-attribute {
	color: #eab700;
}

/* Tomorrow Green */
.hljs-string,
.hljs-symbol,
.hljs-bullet,
.hljs-addition {
	color: #718c00;
}

/* Tomorrow Blue */
.hljs-title,
.hljs-section {
	color: #4271ae;
}

/* Tomorrow Purple */
.hljs-keyword,
.hljs-selector-tag {
	color: #8959a8;
}

.hljs {
	display: block;
	overflow-x: auto;
	color: #4d4d4c;
	padding: 0.5em;
}

.hljs-emphasis {
	font-style: italic;
}

.hljs-strong {
	font-weight: bold;
}
</style>

<style>
/*
 * Markdown PDF CSS
 */

 body {
	font-family: -apple-system, BlinkMacSystemFont, "Segoe WPC", "Segoe UI", "Ubuntu", "Droid Sans", sans-serif, "Meiryo";
	padding: 0 12px;
}

pre {
	background-color: #f8f8f8;
	border: 1px solid #cccccc;
	border-radius: 3px;
	overflow-x: auto;
	white-space: pre-wrap;
	overflow-wrap: break-word;
}

pre:not(.hljs) {
	padding: 23px;
	line-height: 19px;
}

blockquote {
	background: rgba(127, 127, 127, 0.1);
	border-color: rgba(0, 122, 204, 0.5);
}

.emoji {
	height: 1.4em;
}

code {
	font-size: 14px;
	line-height: 19px;
}

/* for inline code */
:not(pre):not(.hljs) > code {
	color: #C9AE75; /* Change the old color so it seems less like an error */
	font-size: inherit;
}

/* Page Break : use <div class="page"/> to insert page break
-------------------------------------------------------- */
.page {
	page-break-after: always;
}

</style>

<script src="https://unpkg.com/mermaid/dist/mermaid.min.js"></script>
</head>
<body>
  <script>
    mermaid.initialize({
      startOnLoad: true,
      theme: document.body.classList.contains('vscode-dark') || document.body.classList.contains('vscode-high-contrast')
          ? 'dark'
          : 'default'
    });
  </script>
<h1 id="%E6%97%B6%E9%97%B4%E5%BA%8F%E5%88%97%E4%B8%8E%E9%9A%90-markov-%E6%A8%A1%E5%9E%8B">时间序列与隐 Markov 模型</h1>
<p>[TOC]</p>
<p>之前讨论的模型都是在独立样本上进行学习。如果样本不独立，那么学习任务会变成什么样呢？这就是本章要讨论的课题。</p>
<p><strong>约定</strong> 出于书写上的便利和美观，时间序列有时简单记为$X_1 \cdots X_T$。序列$X_{s}\cdots X_{t}={X_{s},X_{s+1},\cdots, X_{t}}$进一步简记为$X_{s:t}$。一般地，若指标集$I={t_1, \cdots, t_M}\subset{1,\cdots,T}$，则$X_I$表示$X_{1:T}$的子序列$X_{t_1}\cdots X_{t_M}$。最后令$X_{I^c}={X_t,t\notin I}$代表$X_I$的补序列。（记号$s:t$的灵感来自计算机语言的数组索引或切片格式。笔者认为这个记法简洁又不失严谨。）</p>
<h2 id="%E6%97%B6%E9%97%B4%E5%BA%8F%E5%88%97%E7%AE%80%E4%BB%8B">时间序列简介</h2>
<p>时间序列是一个随时间变化的随机变量，记为${X_t}$，其中$t=1,2,\cdots$。时间序列在统计学中就已经被充分研究了。因为序列中的随机变量不一定独立，所以对于长度为$T$的序列，要用联合分布来描述：
$$
P(X_1,X_2,\cdots, X_T)
$$
对无限序列，我们必须用一系列联合分布来描述，
$$
P(X_1,X_2,\cdots, X_t), t=1,2,\cdots
$$
因此无限长序列可以归结为任意有限长序列，而我们也就不在乎序列是否有限了。这些联合分布统称为序列分布。</p>
<p><em>例</em> 最典型的时间序列模型是<strong>自回归滑动平均模型(ARMA)</strong>：
$$
X_t = \sum_{i=1}^p a_iX_{t-i} + Z_t + \sum_{i=1}^q a_iZ_{t-i},
$$
其中$Z_i\sim N(0,\sigma^2)$独立同分布，记$X_t\sim ARMA(p,q)$。</p>
<p>下面这个定义只是为本章的论述提供一个根据。</p>
<p><strong>定义（时间序列）</strong>
时间序列是指由服从某个给定分布族$\mathcal{F}$的随机变量构成的序列${X_t,t=1,2,\cdots}$。因此$X_t$有共同的样本空间$\mathcal{X}$，而且$X_t\sim P_{t}\in\mathcal{F}$，但分布$P_t$不一定相同。当所有$X_t$是同分布时，即$X_t\sim P_t=P$，序列${X_t}$也是分布$P(X)$的样本，只是存在相关性。此时，可称$P$为<strong>背景分布</strong>，以和序列分布相区别。（根据符号约定，也可记为${X_t}\sim P^*$。）</p>
<p><em>注</em> 对于时间序列模型，无论所有随机变量$X_t,t=1,\cdots, T$是否同分布的，都不能直接用最大化“独立似然”$\sum_t\log P(X_t)$进行估计参数。</p>
<h3 id="markov-%E9%93%BE">Markov 链</h3>
<p>如果时间序列中的$X_t$和$X_{1:t-h-1}$相对$X_{t-1:-th}$条件独立，则被称为$h$-阶 <strong>Markov链</strong>。本章只关注1-阶Markov链，其确切定义如下。
<strong>定义（Markov链）</strong>
Markov链的联合分布为
$$
P(X_1,\cdots,X_T) = P(X_1)\prod_{t=1}^{T-1} P(X_{t+1}|X_{t})
$$
即$X_{t+1}$在$X_t$的条件下和其他变量无关。当条件分布$P(X_{t+1}|X_{t})$和时间无关时，称该Markov链是<strong>齐次的</strong>。概率图模型为</p>
<pre class="hljs"><code><div>x-&gt;x...-&gt;x
</div></code></pre>
<p><em>注</em> 注意$P(X_{t+1}|X_{t})$是习惯写法，在齐次条件下，作为二元函数，是和时间无关的。</p>
<p>齐次Markov链是我们真正感兴趣的。齐次Markov链的联合分布完全由和时间无关的初始概率$P(X_1)$和转移概率$P(X_{t+1}|X_{t})$决定。我们可以用三元组
$$
(\mathcal{X}, P(X_1),P(X_{t+1}|X_{t}))
$$</p>
<p>来表示这样的序列，其中$P(X_1)$称为<strong>初始概率</strong>，$P(X_{t+1}|X_{t})$称为<strong>转移概率</strong>。序列样本${X_t}$的似然为，
$$
l(p)=\sum_{t=1}^{T-1}\log p(X_{t+1}|X_{t}) +\ln p(X_1).
$$
若观测到多个独立序列${X_{i,t}, t=1\cdots,T_i},i=1,\cdots N$，则相应的似然为，
$$
l(p)=\sum_i(\sum_{t=1}^{T_i-1}\log p(X_{i,t+1}|X_{i,t}) +\ln p(X_{i,1})).
$$</p>
<p><em>注</em> 从序列分布的角度讲，一个样本序列${X_t}$只能算一个样本点。若干个序列才构成一个序列分布的样本。</p>
<p><strong>事实</strong>
齐次Markov链的参数估计等价于判别模型$P(X_{t+1}|X_t)$基于独立样本${(X_t,X_{t+1}),t=1,\cdots, T-1}$的监督学习，带先验项$P(X_1)$。</p>
<p>出现在Markov链的样本相关性也被称为<strong>时间依赖性</strong>。下一节，我们将把样本的时间依赖性装配到机器学习中。</p>
<p>一般的确定型动力系统由状态映射$x_{t+1}=f(x_t,t)$和初始状态$x_1$表示。Markov 链显然是这类系统的概率模型。Markov链有时也给出初始状态。此时考虑条件似然
$$
l(p;X_{2:T}|X_1)=\sum_{t}\log p(X_{t+1}|X_{t}).
$$</p>
<p>根据上述讨论，我们构造一个基于时间序列的决策模型。
<strong>Markov决策模型</strong>
Markov决策模型由下述经验风险确定，
$$
J(\theta):=\sum_{t=1}^{T-1} w_tl(\theta;X_{t},X_{t+1}) + l(X_1,\theta)
$$
其中$w_t$是合理的时间权重，${X_t}$是样本序列。</p>
<h3 id="%E7%8A%B6%E6%80%81%E7%BC%BA%E5%A4%B1-markov-%E9%93%BE">状态缺失 Markov 链</h3>
<p>我们来考虑一个简单的隐变量模型。
<strong>定义（状态缺失 Markov链）</strong>
状态缺失 Markov链有和Markov链一样的联合分布，其中只观测到样本子序列$X_I$的值。</p>
<p>设$I={i_1,\cdots,i_m}$。于是状态缺失 Markov链的边缘分布为：
$$
p(X_I)=\sum_{t_i,t_{i+1}\in I}\log p(X_{t_{i+1}}|X_{t_i}) +\ln p(X_{i_1}).
$$
因此，状态缺失 Markov链（边缘分布确定的Markov链）依然是Markov链；而状态缺失齐次 Markov链可能不是齐次Markov链。</p>
<h2 id="%E6%97%B6%E9%97%B4%E5%BA%8F%E5%88%97%E7%9A%84%E7%9B%91%E7%9D%A3%E5%AD%A6%E4%B9%A0">时间序列的监督学习</h2>
<p>现在设想从生成模型$p(X,Y)$中采样：先得到样本${y_t}$；再独立地在$y_t$条件下生成$x_t$；最终得到样本$(x_t,y_t),t=1,2,\cdots$。记输入和输出空间记为$\mathcal{X}, \mathcal{Y}$。若样本${y_t}$是时间依赖的(齐次 Markov链)，则我们构造了一个关于时间序列${(x_t,y_t),t=1,2,\cdots}$的模型，其转移概率为（与时间$t$无关），
$$
p(x_{t+1},y_{t+1}|x_t, y_{t})=p(y_{t+1}| y_{t})p(x_{t+1}|y_{t+1})
$$
​​​起始概率为，$p(y_1)p(x_1|y_1)$。该序列模型的严格定义如下。</p>
<p><strong>定义（作为机器学习的 Markov 模型-生产形式）</strong>
齐次Markov时间序列监督学习生成模型（简称 Markov 模型）可表示为一个五元组
$$
(\mathcal{X}, \mathcal{Y}, p(y_1), p(x_t|y_t), p(y_{t+1}|y_t)),
$$
其中$p(y_1), p(y_{t+1}|y_t), p(x_t|y_t)$分别是模型的起始概率，转移概率和<strong>发射概率</strong>（即生成概率）。样本的依赖关系如前所述，另见图X.X。</p>
<p><img src="../src/TS.png" alt=""></p>
<p><em>注</em> 再次强调，发射概率和转移概率和时间$t$无关，形式上不含$t$，换言之$P(x_t|y_t)$和$P(y_{t+1}|y_{t})$各自是同分布的（在$y_t=y_t$的条件下）。</p>
<p>对于 Markov 模型，我们只能考虑样本序列${(x_t,y_t)}$上的似然：
$$
l(p)=\ln p(x_{1:T},y_{1:T})\
={\ln p(y_1)+\sum_{t=1}^{T-1}\ln p(y_{t+1}|y_{t})+\sum_{t=1}^T\ln p(x_t|y_t)}
$$</p>
<p>(1)包含三个似然函数：初始似然$\ln p(y_1)$, 生成（发射）似然$\sum_{t=1}^T\ln p(x_t|y_t)$，还有转移似然$\sum_{t=1}^{T-1}\ln P(y_{t+1}|y_{t})$。</p>
<p>正如普通统计学习模型是一种特殊的、有内在结构的统计模型，Markov 模型是一种特殊的，有内在结构的齐次Markov链。</p>
<p><strong>事实</strong>
定义X.X 的 Markov 模型等价于转移概率为$p(x_{t+1},y_{t+1}|x_t, y_{t})=p(y_{t+1}| y_{t})p(x_{t+1}|y_{t+1})$的齐次Markov链${(x_t, y_{t}),t=1,\cdots,T}$。</p>
<p>从似然(X.X)看，Markov 模型是序列到序列的监督学习模型。换个角度看（如本节开头所述），它也是监督学习模型构成的序列。</p>
<p><strong>事实</strong>
Markov 模型的参数估计等价于生成模型$p(x_t,y_t)$基于独立样本${(x_t,y_t),t=1,\cdots, T}$的监督学习，带正则项$p(y_1)\prod_t p(y_{t+1}|y_{t})$，也等价于判别模型
$$
p(y_{t+1}|y_t)p(x_{t+1}|y_{t+1})
$$
基于独立样本${(y_t,x_{t+1},y_{t+1}),t=1,\cdots, T-1}$的监督学习，带正则项$p(y_1)p(x_1|y_1)$，其中$(y_t,x_{t+1})$是输入，$y_{t+1}$是输出。</p>
<p><em>注</em> Markov 模型的参数估计也等价于生成模型$p(x_t,y_t)$基于独立样本${(x_t,y_t),t=1,\cdots, T}$的监督学习，带正则项$p(y_1)\prod_t p(y_{t+1}|y_{t})$。不过，这个正则项看上去有点繁琐。</p>
<p><em>注</em> 状态自动机由条件状态映射$y_{t+1}=f(y_t,x_{t+1})$和初始状态$y_0$表示。Markov 模型显然是这类系统的概率模型。</p>
<p>我们也可以定义Markov模型的判别形式，并将关于它的讨论留给读者。
<strong>定义（Markov 模型-判别形式）</strong>
Markov 模型判别形式或 Markov 判别模型的联合分布为：
$$
p(y_{1:T}|x_{1:T})= \prod_{t=1}^{T}p(y_t|x_t) \prod_{t=1}^{T-1} p(y_{t+1}|y_t)
$$
其中${x_t}$是给定的（不妨称为“设计序列”）。该模型可表示为一个四元组
$$
(\mathcal{X}, \mathcal{Y}, p(y_t|x_t), p(y_{t+1}|y_t)),
$$
其中$p(y_{t+1}|y_t), p(y_t|x_t)$分别是模型的转移概率和<strong>判别概率</strong>。</p>
<p><em>注</em> Markov模型的另一种定义是：没有$y_1$所对应的$x_1$（或$x_1$不可观测），或者在$y_1$前加入$y_0$（没有对应的$x_0$）。这么做会使联合分布的表达式更。</p>
<p><em>注</em> 无论多复杂的序列模型，都是样本序列空间$(\mathcal{X}\times\mathcal{Y})^*$上的统计模型。</p>
<p><strong>约定</strong> 遵守约定X.X。若观测到多序列样本${(x_{i\cdot}, y_{i\cdot}), x_{i\cdot}=x_{i1}\cdots x_{iT_i},y_{i\cdot}=y_{i1}\cdots y_{iT_i}, i=1,\cdots,N}$，包含$N$个序列，且第$i$个序列长度为$T_i$，则模型的似然为，
$$
l(p)=\sum_i\ln p(x_{i\cdot},y_{i\cdot})\
=\cdots
$$</p>
<h3 id="%E6%A0%87%E5%87%86-markov%E6%A8%A1%E5%9E%8B">标准 Markov模型</h3>
<p>现在只考虑Markov模型的参数形式：$p(x_{1:T},y_{1:T}|\theta)$。</p>
<p><strong>定义 Markov序列模型的参数解绑形式（基础型）</strong>
将Markov序列模型的参数$\theta$分解成三部分$\theta=(\pi,\alpha,\beta)$，每部分独立表征不同的似然，即
$$
l(\pi,\alpha,\beta)=\ln p(y_1|\pi)+\sum_{t=1}^{T-1}\ln p(y_{t+1}|y_{t},\alpha)+\sum_{t=1}^T\ln p(x_t|y_t,\beta).
$$</p>
<p>显然，参数解绑情况下，联合似然归结为三个有独立参数的似然$\ln p(y_1|\pi)$, $\sum_{t=1}^{T-1}\ln p(y_{t+1}|y_{t},\alpha)$, $\sum_{t=1}^T\ln p(x_t|y_t,\beta)$，其中前两个是齐次Markov链的似然，而第三个是生成模型$p(x_t|y_t)$关于独立样本的生成似然（独立生成似然）。</p>
<p><strong>事实</strong>
参数解绑形式(X.X)使Markov序列模型退化为两个独立的模型：一个生成模型和一个Markov链，其中生成模型不含先验分布。</p>
<p>注X.X 已经指出一个样本序列只能算一个样本点。若观测到多个独立样本序列${(x_{i\cdot}, y_{i\cdot}), x_{i\cdot}=x_{i1}\cdots x_{iT_i},y_{i\cdot}=y_{i1}\cdots y_{iT_i}, i=1,\cdots,N}$，则三部分参数的MLE等价于独立地优化下述三个似然（多序列形式）：
$$
\hat{\pi}=\argmax_{\pi} \sum_i\ln p(y_{i1},\pi)\
\hat{\alpha}=\argmax_{\alpha}\sum_i\sum_{t=1}^{T_i-1}\ln p(y_{i,t+1}|y_{it},\alpha)\
\hat{\beta}=\argmax_{\beta}\sum_{k}\sum_{y_{it}=k}\ln P(x_{it}|k,\beta)
$$</p>
<p>从现在开始我们只关心$y_t$是离散的情形（取$K$个值）。
<strong>定义（标准 Markov模型）</strong>
标准 Markov模型可用下述五元组表示：
$$
(\mathcal{X},\mathcal{Y}, \pi, A, \beta)
$$
其中$A_{kl}=p(y_{t+1}=l|y_t=k)$是离散随机变量$y_t$的转移矩阵(大小为$K\times K$)​​，$K$维向量$\pi$​代表起始分布，即$P(y_1=k|\pi)=\pi_{k}$。</p>
<p>若无特别说明，本书中的“Markov模型”专指这种情形。Markov模型似然函数形式如下：
$$
l(\pi,A,\beta)={\ln \pi_{y_1}+\sum_{t=1}^{T-1}\ln A_{y_{t}y_{t+1}}+\sum_{t=1}^T\ln p(x_t|y_t,\beta)},
$$</p>
<p>对观测到的多序列样本${(x_{i\cdot}, y_{i\cdot}), i=1,\cdots,N}$有如下参数估计：</p>
<p>$$
\hat{\pi}<em>k\sim N(y</em>{i1}=k)\
\hat{A}<em>{lk}\sim N(y</em>{it}=k,y_{i,t+1}=l)\
\hat{\beta}=\argmax_{\beta}\sum_{k}\sum_{y_{it}=k}\ln p(x_{it}|k,\beta)
$$</p>
<p>在全面的参数解绑形式中，参数$\beta$可进一步分离，即$p(\cdot|k,\beta)=p(\cdot|\beta_k)$，那么有
$$
\hat{\beta}<em>k=\argmax</em>{\beta_k}\sum_{y_{it}=k}\ln p(x_{it}|\beta_k)
$$</p>
<h3 id="%E7%89%B9%E6%AE%8A-markov%E6%A8%A1%E5%9E%8B">特殊 Markov模型</h3>
<p><strong>定义 基于特定分布族的 Markov模型</strong>
一般地，如果某个Markov模型的发射概率$P(x_t|y_t)$均服从某个分布族$\mathcal{F}$，那么就称该Markov 模型为$\mathcal{F}$-Markov 模型。</p>
<p><strong>定义 Categorial Markov模型</strong>
设$x_t$也是离散的，服从Categorial 分布(取$L$个值)。令发射概率$p(x_t|y_t)=B_{y_tx_t}$​​，其中$B$是（大小为$K\times L$）随机矩阵​。此模型可用五元组$(\mathcal{X},\mathcal{Y},\pi,A,B)$表示，且有似然函数，
$$
l(\pi,A,B)={\ln \pi_{y_1}+\sum_{t=1}^{T-1}\ln A_{y_{t}y_{t+1}}+\sum_{t=1}^T\ln B_{x_ty_t}}
$$</p>
<p>容易给出 Categorial Markov模型中参数$B$的估计：
$$\hat{B}<em>{kl}\sim N(x</em>{it}=l,y_{it}=k)
$$</p>
<p><strong>定义（Gaussian Markov模型）</strong>
Gaussian Markov模型的发射概率为，
$$
P(x_t|y_t=k)\sim N(\mu_k,\Sigma_k),
$$
其中协方差矩阵$\Sigma_k$可以和GMM一样分成四种类型。此模型可用五元组$(\mathcal{X},\mathcal{Z}, \pi, A, {\mu_k,\Sigma_k})$表示。</p>
<p>Gaussian Markov模型的参数$\mu_k,\Sigma_k$的估计如下：
$$
\hat{\mu}<em>k=\frac{1}{N(y</em>{it}=k)}\sum_{y_{it}=k}x_{it}\
\hat{\Sigma}<em>k=\frac{1}{N(y</em>{it}=k)-1}\sum_{y_{it}=k}(x_{it}-\hat{\mu}<em>k)^T(x</em>{it}-\hat{\mu}_k)
$$</p>
<h3 id="%E5%8F%91%E5%B0%84%E7%BC%BA%E5%A4%B1markov%E6%A8%A1%E5%9E%8B">发射缺失Markov模型</h3>
<p>设脚标集$I\subset {1,\cdots,N}$。若Markov模型的输入子序列$x_{I^c}$没有被观测到，或者说$y_{I^c}$没有哦发射$x_{I^c}$，则所有可观测变量的边缘分布为，
$$
p(x_I, y_{1:T})
\sim p(y_1)\prod_tp(y_{t+1}|y_t)\prod_{t\in I}p(x_t|y_t)
$$</p>
<p>虽然从模型的设定看，这似乎是一个隐变量模型，但是边缘分布(X.X)有一个闭合形式，且只会比Markov模型的联合分布简单。于是引入下述模型。</p>
<p><strong>定义 发射/输入/数据缺失Markov模型</strong>
模型有分布(X.X)刻画，其中参数依然是三个基本概率。为了完整描述该模型，需要加上脚标集$I$（称为“发射脚标集”）。可构造相应的参数分离模型。</p>
<p><em>注</em> 更严谨的分布表达式是$p(x_{1:T}, y_{1:T}|I)=p(x_I, y_{1:T})$，其中发射脚标集$I$是条件变量。</p>
<p>以发射缺失 Categorial Markov模型为例，参数$B$的估计为，
$$\hat{B}<em>{kl}\sim N(x</em>{it}=l,y_{it}=k,t\in I_i)
$$
其中第$i$个序列样本的发射脚标集为$I_i$。</p>
<p>发射缺失Markov模型并不实用，因为它描述的不是常见的情形：在没有观测到$x_t$的情况下，还能观测到$y_t$。因此，我们不会深入讨论这种模型。</p>
<h3 id="%E9%A2%84%E6%B5%8B%E4%B8%8E%E7%94%9F%E6%88%90">预测与生成</h3>
<p>对于参数解绑模型来说，学习过程是可以归结为无时间依赖的机器学习问题，因此不会构成我们的障碍。然而，它的预测或生成任务却相当困难。这使得我们不得不花一定篇幅讲清楚任务的具体内容。由于历史原因，我们把如何完成任务留到下一节讲解。</p>
<h4 id="%E9%A2%84%E6%B5%8B%E4%BB%BB%E5%8A%A1">预测任务</h4>
<p>考虑预测任务时，我们总是假定联合分布$p(x_{1:T},y_{1:T})$是已知的。笔者就懒得用参数形式了。</p>
<p>预测任务的一般形式是，
$$\max_{y_{1:T}}p(y_{1:T}|C)$$，
其中$C$代表某种条件。</p>
<p>先来看两类基本的Markov模型预测任务：</p>
<ol>
<li>$\max_{y_{1:T} }p(y_{1:T}|x_{1:T})$​</li>
<li>$\max_{y_{{T+1}:{T+S}}} p(y_{{T+1}:{T+S}}|y_{1:T})$或 $\max_{y_{{T+1}:{T+S}}} p(y_{{T+1}:{T+S}}|x_{1:T})$​</li>
</ol>
<p>预测任务1是普通统计学习的预测或判别任务；预测任务2是时间序列模型特有的时间预测任务。此外，预测任务1似乎要在$K^T$大小的集合中所有最优解。当$T$很大时（通常如此），这是不可承受的。因此我们要构造高效的算法。</p>
<p><strong>定义 数据缺失条件下的预测任务/数据缺失问题</strong>
若输入子序列$x_{I^c}$没有被观测到，则只能用子序列$x_I$预测$z_{1:T}$。该任务归结为求解下述优化问题，
$$
\max_{y_{J^c}}{p(y_{1:T}|x_I)
\sim p(y_1)\prod_tp(y_{t+1}|y_t)\prod_{t\in I}p(x_t|y_t)
}
$$</p>
<p>数据缺失条件下的预测任务就是发射缺失Markov模型的预测。显然，这样的任务在普通机器学习模型中是不成立的，因为样本是独立，自然不能预测$z_t,t\in I^c$。</p>
<p>在部分输出$y_t,t\in J$是可以观测到的情况下，预没有被观测到的$y_{t},t\notin J$。具体地说，我们要求解下述优化问题，
$$
\max_{y_{J^c}}{p(y_{J^c}|x_I,y_J)\
\sim p(y_1)\prod_{t\in J^c \text{或}t+1\in J^c}p(y_{t+1}|y_t)\prod_{t\in I}p(x_t|y_t)}
$$</p>
<h4 id="%E7%94%9F%E6%88%90%E4%BB%BB%E5%8A%A1">生成任务</h4>
<p>时间序列模型的生成任务一般形式为，
$$
x_{1:T}\sim p(x_{1:T}|C)
$$
其中$C$代表某种条件。根据Markov模型的定义，在得到$y_t$后，可以根据发射概率$p(x_t|y_t)$独立地生成$x_t$。因此，如果没有$y_t$，那么我们总是事先根据响应/预测概率$p(y_{1:T}|C)$生成$y_{1:T}$。</p>
<p>例如，从某个确定的状态$y_1$开始生成$x_{1:T}$​​，即从分布$p(x_{1:T}|y_1)$​中采样。具体见下述算法。</p>
<p><em>基本的生成算法</em>
输入: $y_1$（或从初始分布$p(y_1)$中采样）
返回：$x_{1:T}$</p>
<ul>
<li>$y_2\cdots y_T\sim p(y_1,y_2,\cdots, y_T)$（可以是随机采样，也可以是取众数，也可以事先给定）;</li>
<li>采样$x_t\sim p(x_t|y_t), t=1,\cdots, T$;（也可取众数）</li>
</ul>
<p>上述生成算法指的是生成序列以$y_1$为条件，其中$y_1$是给定的或从初始分布$p(y_1)$中采样。当然算法也可以以$x_1$为条件。更贴合实际地，从$x_{1:T}$​生成$x_{{T+1}:{T+S}}$​​。此时应该先根据预测概率$p(y_T|x_{1:T})$预测$y_T$，再执行生成算法。这就是是“AI续写”的原理，只是Markov模型实现的效果会很差。我们不妨称之为“<strong>续写任务</strong>”。相应的生成算法就称为续写算法。</p>
<p><em>续写算法</em>
输入：样本序列$x_{1:T}$
返回：样本序列$x_{{T+1}:{T+S}}$</p>
<ul>
<li>解$\max_{y_T}p(y_{T}|x_{1:T})$预测$y_T$，或从$p(y_{T}|x_{1:T})$随机采样;</li>
<li>$y_{T+1}\sim p(y_{T+1}| y_{T})$（或取众数）;</li>
<li>根据生成算法X.X，从$y_{T+1}$生成序列$x_{{T+1}:{T+S}}$。</li>
</ul>
<p>笔者特别引入一种生成任务。
<strong>定义 Markov模型的完形任务</strong>：只用输入子序列$x_I$预测整个$x_{1:T}$，其中$I={t_1,\cdots,t_M},1\leq t_1 &lt;\cdots &lt; t_M \leq T$，即在$t_i$和$t_{i+1}$之间填补$x_t,t_i&lt;t&lt;t_{i+1}$。任务的关键是根据下述预测概率采样：
$$
p(y_{1:T}|x_I)\sim p(y_1)\prod_tp(y_{t+1}|y_t)\prod_{t\in I}p(x_t|y_t)
$$</p>
<p>完形任务和数据缺失条件下的预测任务的设定是一样的，但最终的目标不同。完形任务可以利用条件预测任务的结果。完形任务属于生成任务；续写任务是它的特殊形式。
更一般地，考虑输出子序列$y_J$被观测到的完形任务：
$$
y_{J^c}\sim p(y_1)\prod_{t+1\in J^c \text{或}t\in J^c}p(y_{t+1}|y_t)\prod_{t\in I}p(x_t|y_t)}
$$</p>
<p><em>完形算法</em></p>
<ul>
<li>根据预测概率(X.X)或(X.X)采样，得到$y_{1:T}$ (或取众数，即条件预测)</li>
<li>然后根据发射概率$p(x_t|y_t)$逐个生成$x_t,t\notin I$（或取众数）。</li>
</ul>
<h2 id="%E9%9A%90markov%E6%A8%A1%E5%9E%8B">隐Markov模型</h2>
<p>下面我们学习本章“主角”——隐Markov模型（HMM）。HMM广泛应用于序列形式信息的处理，比如语音识别，文本序列，生物基因序列等。在耐心学习它的原理之后，我们会给出几个有趣的应用。</p>
<p>就其结构和算法而言，HMM是本书最复杂的模型之一。需要初学者花更多的时间学习。不过在了解时间序列和隐变量模型之后，学习HMM就不那么困难。简单地说，HMM是Markov 序列模型的无监督版本，属于隐变量模型的范畴。它有和 Markov 序列模型一样的模型和概率图，其中$Z_t$是不可观测的。</p>
<pre class="hljs"><code><div>X_1      ... X_t      X_{t+1}      ...
^            ^        ^
|        ... |        |            ...  
|            |        |            
Z_1 ---&gt; ... Z_t ---&gt; Z_{t+1} ---&gt; ...  
</div></code></pre>
<p>此外，根据 Markov 序列模型的似然函数，可以立刻构造出HMM的Q函数：
$$
Q(p,q)=\int q(z_{1:T}){\ln p(z_1)+\sum_{t=1}^{T-1}\ln p(z_{t+1}|z_{t})+\sum_{t=1}^T\ln p(x_t|z_t)}d{z_{1:T}}\
=\int q(z_1)\ln p(z_1)d{z_1} + \sum_{t=1}^{T-1}\int q(z_t,z_{t+1})\ln p(z_{t+1}|z_{t})dz_t dz_{t+1}\
+\sum_{t=1}^T\int q(z_t)\ln p(x_t|z_t)d{z_t}
$$</p>
<p><em>注</em> 特别强调，边缘分布$q(z_t)$和联合分布$q(z_t,z_{t+1})$可能和$t$是有关的，不仅仅和$z_t$的取值有关。当然，我们可以进一步限制分布$q(z_1,\cdots,z_T)$的形式。</p>
<p>不难发现，Q函数的第三个求和项是重构误差（的相反数）。显然Q函数还可以写成
$$
Q(p,q)=\sum_tQ_t(p,q)
$$
其中
$$
Q_1(p,q):=\int q(z_1)\ln p(x_1,z_1)d {z_1}\
Q_t(p,q):=\int q(z_{t-1},z_{t})\ln p(x_{t},z_{t})p(z_{t}|z_{t-1})d z_{t-1}d z_{t},t=2,\cdots,T
$$</p>
<p>正如在Markov模型中，似然函数可以分解成相对独立的几个部分，在这里，Q函数可以分解成相对独立的几个部分（都是Q函数）。但这个公式可能没有什么用处。</p>
<p>要充分利用变分技巧，肯定要计算F函数。如果变分分布$q$没有特别性质，那么只能按照关系式(X.X)计算。如果$q$也有Markov性，那么可以写出与(X.X)类似形式的F函数：
$$
F(p,q)=\int q(z_1)\ln \frac{p(z_1)}{q(z_1)}d{z_1} + \sum_{t=1}^{T-1}\int q(z_t,z_{t+1})\ln \frac{p(z_{t+1}|z_{t})}{q(z_{t+1}|z_t)}dz_t dz_{t+1}\
+\sum_{t=1}^T\int q(z_t)\ln p(x_t|z_t)d{z_t}
$$
且有(X.X)那样的分解式。</p>
<h3 id="hmm-%E7%9A%84%E6%A0%87%E5%87%86%E5%BD%A2%E5%BC%8F">HMM 的标准形式</h3>
<p>在HMM的语境下，作为隐变量的$Z_t$也称为<strong>状态</strong>。所有可能的状态构成<strong>状态空间</strong>，记为$\mathcal{Z}$​。输入$X_t$的样本空间依然记为$\mathcal{X}$。同样我们也只关注$Z_t$是离散的情形（离散状态）。</p>
<p>HMM 和 Markov 模型的联合分布有一样的形式。我们简要叙述其定义。
<strong>定义 HMM</strong>
HMM由下述五元组刻画。
$$
(\mathcal{X},\mathcal{Z}, \pi, A, \beta)
$$
其中$\mathcal{Z}={1,\cdots,K}$是离散状态空间，也是模型的隐空间，$\pi, A, \beta$分别代表转移概率，发射概率和初始概率。若无歧义，只写出参数$\theta=(\pi, A, \beta)$。</p>
<h4 id="%E5%87%86%E5%A4%87">准备</h4>
<p>除了基本概率（转移概率等），为了设计高效的算法，还要定义两个概率
$$
\gamma_t(k):=p(z_t=k|x_{1:T}, \theta)\
\xi_t(k,l):=p(z_t=k,z_{t+1}=l|x_{1:T},\theta)
$$
其中$t=1,\cdots, T, k,l=1,\cdots,K$。</p>
<p><em>注</em> 两个概率可分别称为$t$时刻响应概率，和$t$时刻转移联合概率，并共同构成模型的充分统计量。作为条件变量的$x$在表达式中总是尽可能被省略。</p>
<p>根据模型定义，HMM 的Q函数为（取变分分布为模型响应分布），
$$
Q(\theta,\theta')
=\sum_{z_1}p(z_1|x_{1:T},\theta')\ln\pi_{z_1}+ \sum_{t=1}^{T-1}\sum_{z_t,z_{t+1}}p(z_t,z_{t+1}|x_{1:T},\theta')\ln A_{z_tz_{t+1}}\
+\sum_{t=1}^T\sum_{z_t}p(z_t|x_{1:T},\theta')\ln p(x_t|z_t,\beta)\
=\sum_{k}\gamma_1(k)\ln\pi_{k}+ \sum_{kl}\sum_{t=1}^{T-1}\xi_t(k,l)\ln A_{kl}\
+\sum_{k}\sum_{t=1}^T\gamma_t(k)\ln p(x_t|k,\beta)
$$</p>
<p>现在可以设想HMM的EM算法大致内容了。不难发现，$A$和$\pi$的更新规则为：
$$
A_{kl} \sim \sum_t\xi_t(k,l),\
\pi_k\sim\gamma_1(k).
$$
不过计算$\gamma_t(k),\xi_t(k,l)$不是件容易的事情。我们将在下一节讨论这一点。此外，$\beta$的估计归结成了无时间依赖的MWLE：
$$
\max_\beta\sum_{kt=1}\gamma_t(k)\ln p(x_t|k,\beta)
$$</p>
<!-- 等价地说，根据独立样本$\{(x_t,k),t=1,\cdots,T, k=1,\cdots,K\}$，最小化下述损函数，
$$
loss(\beta)=-\sum_{kt=1}\gamma_t(k)\ln p(x_t|k,\beta)
$$ 
进一步考虑更彻底的参数解绑形式：$p(x_t|k,\beta)=p(x_t|\beta_k),\beta=\{\beta_k\}$。则$\beta$的估计归结为$K$个独立的MWLE：
$$
\max_{\beta_k}\sum_{t=1}\gamma_t(k)\ln p(x_t|\beta_k)
$$ -->
<h4 id="%E7%89%B9%E6%AE%8A-hmm">特殊 HMM</h4>
<p><strong>定义 基于特定分布族的 HMM</strong>
类似于Markov模型，如果一个HMM的发射概率$P(x_t|z_t)$均服从某个分布族$\mathcal{F}$，那么就称该HMM为$\mathcal{F}$-HMM。</p>
<p><strong>categorical HMM</strong>和categorical Markov模型有相同的联合分布，都用五元组$(\mathcal{Z}, \mathcal{X}, \pi, A, B)$表示，其中$\mathcal{Z}$是隐空间。</p>
<p>利用(4)，立刻写出categorical HMM的Q函数：
$$
Q(\theta,\theta')
= \sum_{k}\gamma_1(k)\ln\pi_{k}+ \sum_{t=1}^{T-1}\sum_{kl}\xi_t(k,l)\ln A_{kl}\
+\sum_{kl}\sum_{x_t=l}\gamma_t(k)\ln B_{kl}
$$
其中$\sum_{x_t=l}$表示对所有满足${x_t=l}$的样本点求和。</p>
<p>类似地，<strong>Gaussian HMM</strong>用五元组$(\mathcal{Z}, \mathcal{X},\pi, A, {\mu_k,\Sigma_k})$，并立刻写出Gaussian HMM的Q函数：
$$
Q(\theta,\theta')
= \sum_{k}\gamma_1(k)\ln\pi_{k}+ \sum_{t=1}^{T-1}\sum_{kl}\xi_t(k,l)\ln A_{kl}\
-\frac{1}{2}\sum_{kl}\sum_{x_t=l}\gamma_t(k)(\ln |x-\mu_k|^2_{\Sigma^{-1}}+\ln \det \Sigma_k)
$$</p>
<p><em>注</em> 由于HMM默认是离散状态的，离散可以专门用来形容$x_t$。因此离散HMM一般指categorical HMM。</p>
<h4 id="hmm%E7%AE%97%E6%B3%95">HMM算法</h4>
<p>和Markov 模型一样，HMM要完成三类任务：学习、预测和生成。我们知道HMM是隐变量模型，因此可以预料到它有比Markov 模型更复杂的学习过程。</p>
<p>在介绍算法前，我们引入一些复杂的概率。这些概率由基本概率递归计算得到。</p>
<p><strong>定义 (前向概率)</strong>
给定观测序列$x_1\cdots x_t$，前向概率定义为，
$$
\alpha_t(k):=p(x_1,\cdots,x_t,z_t=k)
$$</p>
<p>对每个给定$t$，前向概率是关于$z_t=k$的函数。因此​在时间序列中，前向概率这个概念暗含了一个函数序列。</p>
<p>下面用前向概率来计算似然函数。计算的思路是，先递归计算函数序列$\alpha_t(k),t=1,\cdots, T$，然后利用恒等式$p(x_{1:T})=\sum_k\alpha_T(k)$计算似然。</p>
<p><strong>向前算法</strong></p>
<p>输入：参数$\theta$和序列$x_1\cdots x_T$​
返回：$p(x_1\cdots x_T)$​​</p>
<ol>
<li>$\alpha_1(k)=\pi_kp(x_1|z_1=k)$;</li>
<li>$\alpha_{t+1}(k)=(\sum_j\alpha_t(j)A_{jk})p(x_{t+1}|z_{t+1}=k)$​;</li>
<li>$p(x_1\cdots x_T)=\sum_k\alpha_T(k)$​;</li>
</ol>
<p><strong>定义(后向概率)</strong>
给定观测序列$x_1\cdots x_T$，后向概率定义为
$$
\beta_t(k) := p(x_{t+1},\cdots,x_T|z_t=k)
$$</p>
<p>类似于前向概率，后向概率也是一个函数序列。用后向概率计算似然函数的思路和前向概率是类似的：先通过递归计算函数序列，然后用恒等式计算出似然。过程见下。</p>
<p><strong>向后算法</strong></p>
<p>输入：参数$\theta$和序列$x_1\cdots x_T$​​
返回：$p(x_1\cdots x_T)$​​</p>
<ol>
<li>$\beta_T(k)=1$;</li>
<li>$\beta_{t-1}(k)=\sum_l\beta_{t}(l)A_{kl}p(x_{t}|z_t=l)$​;</li>
<li>$p(x_1\cdots x_T)=\sum_k\pi_ip(x_1|z=k)\beta_1(k)$​​;</li>
</ol>
<p><em>注</em> 直观地说，前项概率的计算是从 $X_{1:t}$​前进到$Z_t$; 向后概率的计算则是从 $X_{t:T}$倒退到 $Z_t$​。</p>
<p><em>注</em> 请不要将两个函数和参数$\alpha,\beta$混淆。</p>
<p><strong>事实</strong>
$$
\gamma_t(k)\sim \alpha_t(k)\beta_t(k),\
\xi_t(k, l)\sim \alpha_t(k)\beta_{t+1}(l)A_{kl}p(x_{t+1}|z_{t+1}=j,\theta).
$$</p>
<h4 id="%E5%AD%A6%E4%B9%A0%E7%AE%97%E6%B3%95">学习算法</h4>
<p>Baum-Welch 算法本质上是HMM的Q函数导出的EM算法。有趣的是，该算法在成熟的EM算法提出之前就被提出来了。[]</p>
<p><strong>Baum-Welch 算法</strong></p>
<p>输入: 样本序列$x_1\cdots x_T$​​
返回: 参数$\theta=(\pi,A,\beta)$​​</p>
<ol>
<li>
<p>初始化 $\theta_0 = (\pi_0,A_0,\beta_0)$​</p>
</li>
<li>
<p>循环:</p>
<ul>
<li>
<p>E-步:</p>
<ul>
<li>在$\theta$下利用前向后向算法，计算$\alpha_t(\cdot), \beta_t(\cdot)$;​</li>
<li>利用基本事实，计算$\gamma_t(\cdot), \xi_t(\cdot)$;​</li>
</ul>
</li>
<li>
<p>M-步:
更新$\theta=(\pi,A,\beta)$。首先如下更新$\pi,A$，
$$
\pi_k\sim\gamma_1(k);\
A_{kl} \sim \sum_t\xi_t(k,l);
$$
然后解下述优化问题，作为$\beta$更新值。
$$\max_\beta\sum_{kt}\gamma_t(k)\ln p(x_t|k,\beta)
$$</p>
</li>
</ul>
</li>
<li>
<p>直到收敛</p>
</li>
</ol>
<!-- M-步显然可由从Q函数中导出。算法实现的难点主要集中在$\alpha_t, \beta_t$的计算上。 -->
<p>现在考虑离散情形，也是HMM算法最古老的形式。我们采用Q函数(X.X)。</p>
<p><strong>Baum-Welch 算法（离散情形）</strong></p>
<p>输入: 样本序列$x_1\cdots x_T$​​
返回: 参数$\theta=(\pi,A,B)$​​</p>
<ol>
<li>
<p>初始化 $\theta_0 = (\pi_0,A_0,B_0)$​</p>
</li>
<li>
<p>循环:</p>
<ul>
<li>E-步: 同算法X.X</li>
<li>M-步:
$\pi, A$更新同算法X.X，如下更新$B$，
$$
B_{kl} \sim \sum_{X_t=l}\gamma_t(k);
$$</li>
</ul>
</li>
<li>
<p>直到收敛</p>
</li>
</ol>
<p>M-步显然可由Gibbs不等式从Q函数中导出。算法实现的难点主要集中在$\alpha_t, \beta_t$的计算上。</p>
<h4 id="hmm%E9%A2%84%E6%B5%8B%E4%B8%8Eviterbi-%E7%AE%97%E6%B3%95">HMM预测与Viterbi 算法</h4>
<p>之前就提醒过读者，HMM（或Markov模型）的预测算法并不简单。如常假定参数$\theta$或序列联合分布$p(x_{1:T},z_{1:T})$已知。</p>
<p>这里要利用HMM的联合分布的递归公式：
$$
p(x_{1:T},z_{1:T})=p(x_{1:T-1},z_{1:T-1})p(z_{T}|z_{T-1})p(x_T|z_T)
$$</p>
<p><strong>定义 概率最大值函数</strong>：
$$
\delta_t(k):=\max_{z_{1:t-1}}p(z_{1:t-1},z_t=k,x_{1:T}), \
t=1,\cdots,T,k=1,\cdots,K,
$$
其中$x_{1:T}$为观测到的样本序列。</p>
<p>显然
$$
\max_{z_{1:T}} p(x_{1:T},z_{1:T})=\max_k\delta_T(k)
$$
因此，我们只需计算概率最大值函数$\delta_t(k)$。根据著名的 <strong>Bellman 方程</strong>，可递归计算$\delta_t(k)$，并设计相应的算法。</p>
<p><strong>事实（Bellman 方程）</strong>
$$
\delta_1(k):=p(x_1,z_1=k)\
\delta_{t}(k):=\psi_t(k)p(x_{t} | z_{t}=k)\
\psi_t(k):=\max_l{\delta_{t-1}(l)A_{lk}},t=2,\cdots, T
$$</p>
<p><strong>Viterbi 算法</strong>
输入: $x_1\cdots x_T$
返回: $z_1\cdots z_T$​</p>
<ol>
<li>对每个$k$，​计算$\delta_1(k):=p(x_1,z_1=k)$，</li>
<li>$t$从$2$到$T$执行循环
<ul>
<li>对每个$k$，​计算$\psi_{t-1}(k)$​，</li>
<li>对每个$k$，​计算$\delta_t(k)$，并为每个$k$记录序列$z_1...z_t$，其中$z_t=k$;</li>
</ul>
</li>
<li>最后返回$z_1...z_T, z_T=k$，其中对应$\delta_{T}(k)$最大;</li>
</ol>
<p>为了直观理解 Viterbi 算法，我们构造一个加权有向无环多部图，其中顶点集为${z_t/k, t=1,\cdots, T, k=1,\cdots,K}$（$z_t/k$只是一个形式记号，相当于二元组$(t,k)$），有向边集为$(z_t/k,z_{t+1}/l)$，且该边权重为$\ln p(z_{t+1}=l|z_{t}=k)+\ln p(x_{t+1}|z_{t+1}=l)$。更推荐的做法是，在顶点上分配权重$\ln p(x_{t+1}|z_{t+1}=l)$，而边的权重改为$\ln p(z_{t+1}=l|z_{t}=k)$。见下图。</p>
<!-- [多部图] -->
<p>$\delta_t(k)$对应于从任意第一个子图任一顶点$z_1/\cdot$到第$t$个子图顶点$z_t/k$到最短路径长度（如果顶点有权重，那么路径长度是路径经过的所有边和顶点的权重的和）。在这类图上找最短路径还算比较简单的。</p>
<p>Viterbi 算法是一种动态规划法，可以计算出最优预测，同时避免穷举所有可行解。还有另一种预测算法称为“贪婪算法”。我们在第x章提到决策树就是一种贪婪算法。它通常不能得到最优解，但算法设计简单，计算量少。</p>
<p><strong>贪婪算法</strong>
输入: $x_1\cdots x_T$
返回: $z_1\cdots z_T$</p>
<ol>
<li>计算$z_1=\argmax_{z_1}p(x_1,z_1)$;</li>
<li>计算$z_2=\argmax_{z_2}p(x_2,z_1,z_2)$;</li>
<li>依次计算所有的$z_t$​​​;</li>
</ol>
<h4 id="%E6%95%B0%E6%8D%AE%E7%BC%BA%E5%A4%B1%E4%B8%8B%E7%9A%84%E9%A2%84%E6%B5%8B%E4%B8%8E%E5%AE%8C%E5%BD%A2">数据缺失下的预测与完形</h4>
<p>现在让我们正式完成数据缺失下的预测与完形任务。</p>
<p>回顾数据缺失问题：只观测到$x_t,t\in I$的情况下，预测$z_{1:T}$。</p>
<p><strong>定义 边缘概率最大值函数</strong>：
$$
\delta_t(k):=\max_{z_{1:t-1}}p(x_{I},z_{1:t-1},z_t=k),\
t=1,\cdots,T,k=1,\cdots,K,
$$
其中$x_I:={x_t,t\in I}$为观测到的样本序列。</p>
<p>显然
$$
\max_{z_{1:T}} p(x_I,z_{1:T})=\max_k\delta_T(k)
$$
因此，为了解决数据缺失问题，我们只需计算边缘概率最大值函数$\delta_t(k)$。类似地，构造相应的Bellman 方程和算法。</p>
<p><strong>事实（边缘 Bellman 方程）</strong>
$$
\delta_{1}(k):=\begin{cases}
p(x_{1} | z_{1}=k)\
1,1\not\in I
\end{cases}\
\delta_{t}(k):=\begin{cases}
\psi_t(k)p(x_{t} | z_{t}=k)\
\psi_t(k),t\not\in I
\end{cases}\
\psi_t(k):=\max_l{\delta_{t-1}(l)A_{lk}},t=2,\cdots, T
$$</p>
<p>概括地说，如果$x_{t}$不可观测，那么其发射概率值$p(x_{t} | z_{t}=k)$直接从递归公式中删除。对应的Viterbi 算法做相应调整。用上文中的图论语言说，如果$x_{t}$不可观测，那么对应顶点权重就设置为0。</p>
<p>在预测出$z_{1:T}$后，继续生成$x_{t},t\not\in I$，就能实现完形算法。</p>
<p>最后把状态子序列$z_J$可观测的条件算上。
<strong>定义 条件-边缘概率最大值函数</strong>：
$$
\delta_t(k):=\max_{z_{J^c}}p(x_{I},z_{1:t-1},z_t=k|z_J),\
t=1,\cdots,T,k=1,\cdots,K,
$$
其中$x_I:={x_t,t\in I}$为观测到的样本序列。</p>
<p>相应的 Bellman 方程和Viterbi算法留给读者构造。这里用图论语言描述一下：如果某个$z_t$可被观测到，那么删除第$t$个子图的顶点，只留下$z_t$，且顶点权重设置为0。如果$z_{t+1}$也可以观测到，那么边$(z_t, z_{t+1})$的权重也设置为0。</p>
<h3 id="hmm%E5%8F%82%E6%95%B0%E8%A7%A3%E7%BB%91%E5%BD%A2%E5%BC%8F">HMM参数解绑形式</h3>
<p>在参数解绑形式中，发射概率定设置
$$
P(x_t|z_t=k)\sim p(x|\beta_{k})
$$
模型参数表示为$\theta=(\pi, A, {\beta_k})$​。</p>
<p>于是有
$$
Q(\theta,\theta')
= \sum_{k}\gamma_1(k)\ln\pi_{k}+\sum_{kl}\sum_{t=1}^{T-1}\xi_t(k,l)\ln A_{kl}\
+\sum_{k}\sum_{t=1}^T\gamma_t(k)\ln p(x_t|\beta_k)
$$</p>
<p>参数$\theta$的三部分依然可以按照三个求和项独立更新。此外，$\beta_k$的更新转化为下述若干独立的加权似然的优化问题：
$$
\sum_{t=1}^T\gamma_t(k)\ln p(x_t|\beta_k),k=1,\cdots,K
$$</p>
<p>常见HMM都包含在这个一般形式中。显然离散 HMM符合这种形式, 其中发射概率$p(x_t|z_t=k)\sim Cat(\beta_{k})$​​​。协方差矩阵为 full 类型的 Gaussian HMM 也符合这种形式，其中$p(x_t|z_t=k)\sim N(\mu_{k},\Sigma_{k})$​​。这类 Gaussian HMM的学习算法的M步为(参数${\mu_{k},\Sigma_{k}}$的更新)：
$$
\mu_{k} \leftarrow \frac{\sum_t\gamma_t(k)x_t}{\sum_t\gamma_t(k)}\
\Sigma_{k} \leftarrow \frac{\sum_t\gamma_t(k)(x_t-\mu_k)(x_t-\mu_k)^T}{\sum_t\gamma_t(k)}
$$</p>
<p>这个更新公式和GMM形式上是一样的。</p>
<!-- #### 混合Gauss HMM

混合 Gauss HMM的基本设置为，$x_t|z_t$服从混合Gauss分布。这个模型留给读者研究。 -->
<h3 id="%E5%A4%9A%E5%BA%8F%E5%88%97%E6%83%85%E5%BD%A2">多序列情形</h3>
<p>为了表述简练，上文只考虑单序列情形。把结论推广到多序列情形是容易的：
$$
Q(\theta,\theta')=\sum_{ik}\gamma_1^i(k)\ln\pi_{k}+ \sum_{ikl}\sum_{t=1}^{T_i-1}\xi_t^i(k,l)\ln A_{kl}\
+\sum_{ik}\sum_{t=1}^{T_i}\gamma_t^i(k)\ln p(x_{it}|k,\beta)
$$
其中$x_{i\cdot}={x_{it}}$是观测到的第$i$个序列，长度为$T_i$，相应的辅助概率记为$\xi^i_{t},\gamma^i_t$。</p>
<p>在多序列设定下，Baum-Welch 算法 M-步的更新规则为：
$$
A_{kl} \sim \sum_i\sum_t\xi_t^i(k,l);\
\pi_k\sim\sum_i\gamma_1^i(k);
$$
而$\beta$更新值为下述优化问题的解(MWLE)，
$$
\max_\beta\sum_{ik}\sum_{t=1}^{T_i}\gamma_t^i(k)\ln p(x_{it}|k,\beta)
$$
若参数$\beta$可继续分离，$p(x_t|k,\beta)=p(x_t|\beta_k)$，则$\beta_k$是下述优化问题的解，
$$
\max_{\beta_k}\sum_{i}\sum_{t=1}^{T_i}\gamma_t^i(k)\ln p(x_{it}|\beta_k)
$$</p>
<h3 id="hmm%E7%9A%84mcemsem%E7%AE%97%E6%B3%95">HMM的MCEM/SEM算法</h3>
<p>作为一个隐变量模型，HMM当然可以利用MCEM/SEM算法。这里只写出HMM的SEM算法（单序列形式）。</p>
<p><strong>算法 HMM-SEM算法</strong></p>
<ol>
<li>初始化HMM参数$\theta_0$</li>
<li>
<ul>
<li>预测步：随机采样/随机预测（调用算法X.X）$z_{1:T}\sim p(z_{1:T}|x_{1:T},\theta_t)$;</li>
<li>估计步：更新$\theta_{t+1}=\hat{\theta}$，其中$\hat{\theta}$是参数$\theta$在样本$(x_{1:T},z_{1:T})$下的估计，即基于Markov模型的监督学习;</li>
</ul>
</li>
<li>重复2直到收敛</li>
</ol>
<p><em>注</em> 如果用Viterbi算法直接预测$z_{1:T}$，那么将得到HMM的硬 EM算法。</p>
<p><strong>算法 随机预测</strong></p>
<ul>
<li>$z_{1}\sim p(x_{1}|z_{1})p(z_1)$;</li>
<li>$t$从$2$到$T$执行循环：
设$z_{1:t-1}$已生成，则$z_t\sim p(z_{t}|z_{t-1})p(x_t|z_t)$</li>
</ul>
<p>SEM算法不依赖Q函数的形式，因此也就不再需要计算向前/向后概率。
此外，随机预测又使我们避免执行Viterbi算法。这使得SEM算法能够处理连续状态HMM。但是当$T$很大时，它会消耗很多时间。</p>
<h3 id="%E4%BB%A3%E7%A0%81">代码</h3>
<h4 id="hmmlearnpython">hmmlearn/Python</h4>
<p>hmmlearn是Python提供的HMM实现，完全遵循scikit-learn API设计原则。</p>
<pre class="hljs"><code><div><span class="hljs-meta">&gt;&gt;&gt; </span><span class="hljs-keyword">from</span> hmmlearn.hmm <span class="hljs-keyword">import</span> MultinomialHMM, GaussianHMM
<span class="hljs-meta">&gt;&gt;&gt; </span>remodel = GaussianHMM(n_components=<span class="hljs-number">3</span>, covariance_type=<span class="hljs-string">"full"</span>, n_iter=<span class="hljs-number">100</span>)
<span class="hljs-meta">&gt;&gt;&gt; </span>remodel.fit(X)  <span class="hljs-comment"># X: n_samples * n_features</span>
GaussianHMM(algorithm=<span class="hljs-string">'viterbi'</span>,...
<span class="hljs-meta">&gt;&gt;&gt; </span>Z = remodel.predict(X) <span class="hljs-comment"># n_samples * n_components_</span>

<span class="hljs-meta">&gt;&gt;&gt; </span>MultinomialHMM(n_components=<span class="hljs-number">2</span>)                            
<span class="hljs-comment"># MultinomialHMM(algorithm='viterbi',...</span>
</div></code></pre>
<p>源码解读见 <a href="https://blog.csdn.net/nbu2004/article/details/118992180">hmmlearn source codes</a></p>
<h4 id="%E6%88%91%E7%9A%84%E4%BD%9C%E5%93%81">我的作品</h4>
<p><a href="https://gitee.com/williamzjc/hmmwork">hmm work</a></p>
<h2 id="%E5%85%B6%E4%BB%96%E8%AF%9D%E9%A2%98">其他话题</h2>
<p>这一节，我们灵活地改造一下HMM。希望能给读者更多的启发。</p>
<p>下面，半监督HMM和数据缺失HMM和完形任务的设定是一样的，但是这里参数是未知的，而我们要完成的是学习任务。</p>
<h3 id="%E5%8D%8A%E7%9B%91%E7%9D%A3hmm">半监督HMM</h3>
<p>在第X.X节，我们给出了半监督学习的第一个例子：半监督K-Means算法。在第X.X节，给出了半监督学习基于隐变量模型的一般形式和算法。构造半监督HMM不再是什么棘手的事情。</p>
<p><strong>定义 半监督HMM</strong>
半监督HMM和HMM有相同的模型分布，但是子序列$z_J$是可以观测的。此时模型的隐变量只有补序列$z_{J^c}$。（$J$称为<strong>缺失脚标集</strong>。）</p>
<p><strong>半监督HMM-SEM算法</strong></p>
<ol>
<li>初始化$\theta_0$</li>
<li>
<ul>
<li>预测步：随机采样/随机预测（调用算法X.X）$z_{J^c}\sim p(z_{J^c}|x_{1:T},z_J,\theta_t)$;</li>
<li>估计步：更新$\theta_{t+1}=\hat{\theta}$，其中$\hat{\theta}$是参数$\theta$在样本$(x_{1:T},z_{1:T})$下的估计($z_{1:T}$由$z_J$和$z_{J^c}$合成);</li>
</ul>
</li>
<li>重复2直到收敛</li>
</ol>
<p><strong>算法 随机预测</strong></p>
<ul>
<li>$z_{1}\sim p(x_{1}|z_{1})p(z_1)$ （若$z_1$没有被观测到，即$1\notin J$）;</li>
<li>$t$从$2$到$T$执行循环：
设$z_{1:t-1}$已生成。
<ul>
<li>若$z_t,z_{t+1}$都没有被观测到，即$t,t+1\notin J$，或者$t=T$
则$z_t\sim p(z_{t}|z_{t-1})p(x_t|z_t)$;</li>
<li>若$z_t$都没有被观测到，但$z_{t+1}$观测到，即$t\notin J,t+1\in J$，
则$z_t\sim p(z_{t}|z_{t-1})p(x_t|z_t)p(z_{t+1}|z_t)$;</li>
<li>否则$z_t$取观测值。</li>
</ul>
</li>
</ul>
<h3 id="%E5%8F%91%E5%B0%84%E7%BC%BA%E5%A4%B1hmm">发射缺失HMM</h3>
<p>发射缺失HMM就是发射缺失Markov模型的无监督版本。</p>
<p><strong>定义 发射缺失HMM</strong>
发射缺失HMM和发射缺失Markov模型形式上是一样，但无法观测状态。显然可以进一步定义发射缺失的半监督HMM。</p>
<p>我们直接构造发射缺失HMM的半监督学习算法。</p>
<p><strong>算法 半监督发射缺失HMM-SEM算法</strong>
参考算法X.X。只写出预测步：随机采样/随机预测（调用算法X.X）$z_{J^c}\sim p(z_{J^c}|x_{I},z_J,\theta_t)$; （$I$和$J$分别是发射脚标集和缺失集）</p>
<p><strong>算法 随机预测</strong></p>
<ul>
<li>$z_{1}\sim p(x_{1}|z_{1})p(z_1)$ （$1\notin J$）;</li>
<li>$t$从$2$到$T$执行循环：
设$z_{1:t-1}$已生成。
<ul>
<li>若$t,t+1\notin J$，或者$t=T$，
则$z_t\sim \begin{cases}p(z_{t}|z_{t-1})p(x_t|z_t),t\in I\
p(z_{t}|z_{t-1}),t\notin I\end{cases}$;</li>
<li>若$t\notin J,t+1\in J$，
则$z_t\sim \begin{cases}p(z_{t}|z_{t-1})p(x_t|z_t)p(z_{t+1}|z_t),t\in I\
p(z_{t}|z_{t-1})p(z_{t+1}|z_t),t\notin I\end{cases}$;</li>
<li>否则$z_t$取观测值。</li>
</ul>
</li>
</ul>
<h3 id="%E8%BF%9E%E7%BB%AD%E7%8A%B6%E6%80%81hmm">连续状态HMM</h3>
<p>连续状态HMM的基本困难是，学习和预测都会牵涉大量积分运算。一种解决办法是使用SEM算法。</p>
<h3 id="%E5%8A%A8%E6%80%81%E8%87%AA%E7%BC%96%E7%A0%81%E5%99%A8">动态自编码器</h3>
<p>HMM是一个隐变量模型，因此肯定是一个自编码器。它是样本序列空间$\mathcal{X}^<em>$上的自编码器，但不是样本空间$\mathcal{X}$上的自编码器。这里的关键是，编码器$p(z_t|x_t)$是和时间$t$有关的。如果按照众数$z_t^</em>$确定$x_t$最终的编码结果，那么存在关系式$z_t^*=f(x_t,t)$。</p>
<p><strong>定义 动态自编码器</strong>
我们把编码结果或编码器本身与样本顺序有关的“编码-解码结构”称为<strong>动态自编码器</strong>。</p>
<p>HMM是一种动态自编码器。相对地，称普通自编码器为“静态自编码器”。ELMo/BERT等动态自编码器在文本处理中大放异彩。</p>
<p><em>注</em> 根据注X.X，一个$\mathcal{X}$上的动态自编码器其实也只不过是样本序列空间$\mathcal{X}^*$上的静态自编码器。</p>
<h3 id="%E6%9D%A1%E4%BB%B6hmm">条件HMM</h3>
<p>就结构而言，条件HMM是本书“全景式”模型，综合了所有基础模型的结构。</p>
<p><strong>定义 条件HMM</strong>
条件HMM模型是条件隐变量模型$P(y,z|x)$的序列形式，其联合分布为
$$
P(y_{1:T},z_{1:T}|x_{1:T})=\prod_{t=1}^T P(y_t|z_t)\prod_{t=1}^{T-1}P(z_{t+1}|z_t,x_{t+1})
$$
其中$x_t,y_t,z_t$分别是输入变量（条件变量），输出变量（可观测变量）和隐变量，且样本的时间依赖性由<strong>条件转移概率</strong>$P(z_t|z_{t-1},x_t)$和<strong>条件发射概率</strong>$P(y_{t}|z_{t},x_{t})$决定。</p>
<h2 id="%E6%80%BB%E7%BB%93">总结</h2>
<p>本章介绍了以HMM为主的时间序列模型。引入这些模型的目的是处理文本理解、语音识别等复杂的序列问题。由于要考虑样本的时间依赖性，时间序列模型需要引入迁移概率来描述。HMM不仅考虑时间依赖性，还是一个隐变量模型，因此模型结构和学习算法都变得较为复杂。它的预测算法也不简单，即使基于简单的背景分布。</p>
<p>我们也讨论了作为HMM的有监督版本，Markov模型。这在其他文献里很少被讨论。HMM比Markov模型更符合一般时间序列问题的设定。</p>
<p>Markov模型/HMM可以完成一些非常独特的预测和生成任务。“完形任务”概括了这些任务的一般形式。</p>
<p>我们还构造了半监督HMM和发射缺失HMM。加上完形任务，所有问题概括为，HMM如何在观测到子序列$x_{I}$和$z_J$的条件下的学习和预测。</p>
<p>% 本节给出的HMM算法都是基于标准的EM算法。读者可以自行构造HMM的MCEM/SEM算法、VEM算法。</p>
<hr>
<p><em>练习</em></p>
<ol>
<li>设$x_1,\cdots,x_T$是Markov链，转移概率记为$p(x_{t+1}|x_t)$. 证明恒等式$\prod_{t=1}^{T-1} p(x_{t+1}|x_{t})=p(x_T|x_1)\prod_{t=2}^{T-1} p(x_{t}|x_{t+1},x_1)$.</li>
<li>试说明Markov 判别模型相当于一个由判别模型构成的序列，并且其参数估计归结为某个判别模型基于独立样本的监督学习。Markov 判别模型是否有对应的隐变量模型（部分$y_t$不可观测）？</li>
<li>证明HMM输入序列的联合分布$p(x_1,\cdots,x_T)$可写成如下混合分布的形式：
$$
\sum_{l}\alpha_lf_l(x_1,\cdots,x_T)
$$</li>
<li>验证事实X.X。并讨论，Gaussian HMM 是否相当于 GMM 的齐次 Markov 链，即构造合理的转移概率$p(x_{t+1}|x_t)$，其中$x_t$从混合Gauss分布。</li>
<li>写出混合Gaussian HMM的定义。并设计相应的学习算法。</li>
<li>请推导半监督HMM的Q-函数。</li>
<li>一种特殊的半监督HMM是：部分序列样本可观测到完整的状态序列，而其他序列样本观测完全不到状态序列。这显然是HMM和Markov模型的简单组合。请写出Q-函数，并构造相应的学习算法。</li>
</ol>
<p><em>Ref</em></p>
<p>R. E. Kalman. Mathematical Description of Linear Dynamical Systems, J.S.I.A.M. Control Ser A. 1963, 1(2):152-192</p>
<p>P. J. Brockwell, R. A. Davis. Time Series: Theory and Methods. Springer.</p>
<p>L. Rabiner and B.-H. Juang. Fundamentals of Speech Recognition. Prentice Hall Signal Processing Series, 1993.</p>
<p>L. R. Rabiner and B.-H. Juang, “Speech recognition: Statistical methods,” Encycl. Lang. Linguist., 2006:1–18</p>
<p>L. R. Rabiner. A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE, 1989, 77(2): 257-286.</p>
<p>Daniel Jurafsky &amp; James H. Martin. Speech and Language Processing, 2021.</p>
<p>Gales and Young (2007). “The Application of Hidden Markov
Models in Speech Recognition”, Foundations and Trends in
Signal Processing, 1 (3), 195–304: section 2.2.</p>
<p>Baum, L. E. and T. Petrie. 1966. Statistical inference for
probabilistic functions of finite-state Markov chains. An
nals of Mathematical Statistics, 37(6):1554–1563.</p>
<p>Baum, L. E. and J. A. Eagon. 1967. An inequality with appli
cations to statistical estimation for probabilistic functions
of Markov processes and to a model for ecology. Bulletin
of the American Mathematical Society, 73(3):360–363.</p>
<p>Viterbi, A. J. 1967. Error bounds for convolutional codes and
an asymptotically optimum decoding algorithm. IEEE
Transactions on Information Theory, IT-13(2):260–269.</p>
<p>Inoue,M. and Ueda,N. (2003) Exploitation of unlabeled sequences in Hidden Markov Models. IEEE Trans. Pattern Anal. Mach. Intell., 25, 1570–1581.</p>
<p>Popov, Alexander &amp; Tatyana, Gultyaeva &amp; Uvarov, Vadim. (2016). Training Hidden Markov Models on Incomplete Sequences. 10.1109/APEIE.2016.7806478.</p>

</body>
</html>
